A self-interaction corrected pseudopotential scheme for magnetic and 

strongly-correlated systems. 
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Local-spin-density functional calculations may be affected by severe errors when applied to the 
study of magnetic and strongly-correlated materials. Some of these faults can be traced back to the 
presence of the spurious self-interaction in the density functional. Since the application of a fully 
self-consistent self-interaction correction is highly demanding even for moderately large systems, 
we pursue a strategy of approximating the self-interaction corrected potential with a non-local, 
pseudopotential-like projector, first generated within the isolated atom and then updated during 
the self-consistent cycle in the crystal. This scheme, whose implementation is totally uncomplicated 
and particularly suited for the pseudopotental formalism, dramatically improves the LSDA results 
for a variety of compounds with a minimal increase of computing cost. 
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I. INTRODUCTION 

Among the innumerable successes of the local-density 
(LDA) and local-spin-density (LSDA) approximations 
to density functional theory (DFT) there are also well 
known, systematic failures 1 ' 2 that compromise the accu- 
racy of predictions for a range of properties, especially in 
magnetic and strongly-correlated 3 compounds. 

Typical examples of LSDA failures are the series 
of transition-metal monoxides, 4 improperly described 
in LSDA as either small-gap Mott-Hubbard antiferro- 
magnetic insulators (MnO, NiO) 5,6 or even ferromag- 
netic and non-magnetic metals (FeO, CoO, CuO), 5,7 
whereas according to experiments these materials are 
charge-transfer antiferromagnetic wide-gap insulators. 
A similar situation occurs for the high-T c compounds 
La2CuC>4 and YBa2Cu30g which are described in LSDA 
as nonmagnetic metals instead of as antiferromagnetic 
insulators, 8 and for the perovskite manganites (e.g. 
La^Cai-zMnOs), 9,10 for which the LSDA fails to predict 
the correct magnetic and orbital orderings. In general, 
LSDA favors metallic and ferromagnetic ground states 
over the observed antiferromagnetic insulating ground 
states. This is particularly harmful in the case of hexag- 
onal YMnC>3, which is antiferromagnetic and ferroelec- 
tric, but is described as a metal within LSDA, 11 thus 
preventing the possibility of calculating any ferroelectric 
properties at all. 

These failures can be, at least in part, attributed to the 
presence of the self-interaction (SI) in the LSDA energy 
functional, i.e. the interaction of an electron charge with 
the Coulomb and exchange-correlation potential gener- 
ated by the same electron. The SI vanishes in the thermo- 
dynamic limit for delocalized states, but is present in sys- 
tems characterized by spatially localized electron charges 
like 2p, 3d and 4f electrons. As a consequence of the SI, 
the binding energies, the on-site Coulomb energies (i.e. 
the Hubbard U parameter) and the exchange-splitting 
of the d and f states are underestimated, whereas the 
hybridizations of cation d and anion p states and the 
corresponding band widths are overemphasized. 

Since, for materials with partially filled d states, the 



on-site Coulomb energy should be much larger than the 
charge-transfer energy between p and d electrons, it is 
clear that the suppression of U and the overestimation of 
the p-d hybridization may change dramatically the char- 
acter of the band structure, with a tendency to favor 
metallic and ferromagnetic ground states over insulating 
antiferromagnetic states. Furthermore, within the LSDA 
Kohn-Sham (KS) description based on a local, orbital- 
independent potential, orbital and charge orderings can- 
not be properly accounted for. 

The presence of SI and the possible strategies for elim- 
inating the SI in density functional theories are long- 
standing issues which go back to the Thomas-Fermi and 
the Slater (Xa) approaches. 12 A large amount of lit- 
erature has been produced 13-17 within the context of 
the LSDA. Particularly fundamental are the works by 
Perdew and Zunger 15,17 which proposed a form of self- 
interaction correction (SIC) within LSDA (SIC-LSDA) 
and successfully applied it to the calculation of atomic 
properties. 

More recently, Svane and co-workers presented a long 
series of works 18-28 which attempted the highly challeng- 
ing application of the fully self-consistent SIC-LSDA to 
extended systems with encouraging results. However the 
full SIC-LSDA requires a large computing effort when ap- 
plied to extended systems even for materials with small 
unit cells, and becomes prohibitive for larger systems. 

A useful alternative to this approach has been sug- 
gested by Vogel and co-workers. 29-31 They approximated 
the SI part of the KS potential with a non-local, atomic- 
like contribution included within the pseudopotential 
construction. This scheme, applied to non-magnetic II- 
VI semiconductors and III-V nitrides, is capable of re- 
markable improvements upon the LSDA results while 
keeping the computational cost comparable to that of 
an ordinary pseudopotential calculation. 

Inspired by these results, in this paper we develop an 
approach which is in the same spirit as that of Ref. 29, 
but can be applied to more general cases, and in par- 
ticular to magnetic and highly-correlated systems where 
there is a coexistence of strongly localized and hybridized 
electron charges. The main innovation of our formalism 
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is the introduction of orbital occupation numbers in the 
electron KS Hamiltonian: The pseudopotential-like SIC 
(pseudo-SIC) is rescaled by the occupation numbers cal- 
culated self-consistently within the crystal environment, 
so that the SI coming from localized, hybridized, or com- 
pletely itinerant states may be discriminated accordingly. 
Also, these occupation numbers discriminate between va- 
lence and conduction bands so that only the former are 
corrected, in accord with the idea that conduction states 
are itinerant, thus SI- free. Furthermore, our scheme is 
generally applicable to insulators as well as metals, with- 
out the necessity of knowing the character of the material 
a priori. This is essential for compounds which are er- 
roneously described as metals within LSDA, such as the 
hexagonal YMnC>3, treated in section IV E. 

A second major difference between our approach and 
Ref. 29 lies in the way in which relaxation corrections to 
the SIC potential are accounted for. In Ref. 29 they are 
incorporated as an additional, orbital-dependent number 
which is imported from atomic energy difference calcula- 
tions (i.e. the so-called ASCF method). In contrast, in 
our case the relaxation contribution is included directly 
as an analytical correction in the SIC potential projector, 
and no further atomic calculations are needed. 

In this paper we test the pseudo-SIC on three classes 
of materials, including the wide-band gap insulators (e.g. 
ZnO and GaN), the transition-metal oxides (MnO and 
NiO) and the manganese oxides (YMnOs). Our main fo- 
cus is on the electronic properties (i.e. the band energy 
structure) of these materials, since we want to verify the 
capability of our single-electron Hamiltonian to repro- 
duce the main features of photoemission spectra. For 
some of these compounds we also calculated the theoret- 
ical structure by total energy minimization. In general 
we obtain very encouraging results and clear systematic 
improvements over the LSDA description, without con- 
siderable increase of computing effort. 

Note that we have implemented the pseudo-SIC within 
the ultrasoft pseudopotential method 34 (USPP). This is 
instrumental in keeping a moderate computational cost 
even for large unit cells containing atoms with 2p and 3d 
electrons. 

As is customary when a new scheme is introduced, we 
compare our results with those of other common beyond- 
LDA approaches. According to the set of results pre- 
sented here, the pseudo-SIC ranks among the most ac- 
curate. In particular, it seems to perform equally well 
(or even better) than the very popular LDA+U. 32 Al- 
though a detailed comparison between pseudo-SIC and 
LDA+U is not the aim of this paper, we can point out 
some potential advantages for the pseudo-SIC. First, it 
does not require parameters imported from an external 
theory. Second, it can be applied to both magnetic and 
non-magnetic compounds, whereas the LDA+U is con- 
structed to correct spin-polarized and/or orbital-ordered 
band structures. 

Finally, within LDA+U or even within the exact SIC- 
LSDA a choice of which orbitals arc localized in space and 



thus, which orbitals are to be corrected, has to be made 
before the calculation, and the final result depends on 
this assumption. Instead, in pseudo-SIC the correction is 
applied to all the orbitals indiscriminately, and no choice 
of orbital localization is required. This is an advantage in 
case where a discrimination of the electron charge local- 
ization cannot be established a priori. This happens in 
non-bulk systems such as surfaces and interfaces, when 
localized surface states or resonance states are present, 
or even in bulk materials whenever the charge localiza- 
tion is not, or not only, due to the chemical nature of the 
involved electrons. Examples of this kind are the super- 
conductor cupratcs, in which the d orbital charges have 
different localizations corresponding to different physical 
solutions. 

The remainder of this paper is organized as follows: in 
Section II we review the main features of earlier work on 
SIC implementations in LSDA. In Section III we describe 
the pseudo-SIC formulation applied in combination with 
norm-conserving pscudopotcntials (NCPP). The exten- 
sion of the pseudo-SIC to USPP is given in the Appendix. 
In Section IV we report our results obtained within the 
pseudo-SIC approach, and finally in Section V we present 
the conclusions. 



II. OVERVIEW OF PREVIOUS WORK 

The SI effects in atomic calculations are extensively ex- 
plained in the seminal paper by Perdew and Zungcr. In 
the following we briefly summarize the most fundamen- 
tal features. Due to the SI, the tail of the KS potential 
does not recover its physical long-range limit, -1/r. Thus 
negative ions which are stable experimentally cannot be 
described at all in LSDA since the outermost eigenstates 
are unbound. The atomic total energies are severely over- 
estimated with respect to the experimental values, as a 
consequence of the reduction in binding energy caused 
by the spurious self-screening in the KS potential. 

Consistently, the LSDA KS eigenvalues strongly over- 
estimate the experimental electron removal energies. Al- 
though the KS eigenvalues represent in general Lagrange 
multipliers and the comparison with the removal ener- 
gies may be questioned, calculations 17 clearly show that 
in atoms the SI represents, by far, the largest source of 
mismatch. Furthermore, it can be proved that the high- 
est occupied KS eigenvalue of the exact density functional 
theory equals the atomic ionization potential. 35 Thus, at 
least for this quantity, the discrepancy must be attributed 
to the LSDA inaccuracy. 

The SIC-LSDA proposed by Perdew and Zunger is 
based on the straightforward subtraction of the SI con- 
tribution from the LSDA KS potential: 

VSxefaM - V§ xc [n,m] - V^nf], (1) 

where n and m are the total charge and magnetiza- 
tion densities, and nf is the spin-charge density of the 
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i th orbital. It is understood that in V^ xc [nf] the mag- 
netization density is set equal to 1, i.e. the single elec- 
tron charge is fully spin-polarized. Despite its conceptual 
simplicity, this modification is able to systematically im- 
prove total energies, ionization potentials and electron 
affinities, giving a much better agreement between KS 
eigenvalues and removal or addition energies. 

Notwithstanding this initial success, the SIC-LSDA 
has not enjoyed a large popularity in the first-principles 
community, mainly because the correction severely com- 
promises the feasibility of the LSDA, introducing an 
orbital-dependent, spatially localized (i.e. non-periodic) 
contribution in the KS potential. Thus different wave- 
functions experience different Hamiltonians, and are no 
longer orthogonal. This is not a serious drawback in 
atomic calculations, since the non-orthogonality is a 
small effect and the atomic orbitals can be easily forced 
to be orthogonal during the self-consistent cycle. 17 

In contrast, the direct application of Eq.l to extended 
systems (where the solutions of the KS equations arc 
Bloch states) is very challenging. First of all the SI of 
a Bloch state is an ill-defined quantity that depends on 
the normalization of the wavefunction, and vanishes as 
f2 -1 / 3 in the thermodynamic limit, where O is the sys- 
tem volume. Thus, if we assume that the electron charges 
are properly described as Bloch states the SIC becomes 
immaterial and discardable. 

Clearly this is not necessarily true in general. In cases 
like simple metals or semiconductors with mostly cova- 
lent interactions (e.g. bulk Si) we can safely assume that 
the SI is discardable. However, in many systems the 
electron charges retain atomic-like features such as small 
band dispersion and a spatial distribution strongly lo- 
calized around their ion-cores. A clear example of the 
variable influence of the SI is given by the energy gap 
of semiconductors: for bulk Si the discrepancy between 
the theoretical and experimental energy gaps (0.7 eV and 
1.17 eV, respectively) must be attributed to non-locality 
and many-body effects. Instead, in LiF the SIC accounts 
for almost 95% of the LDA gap error. 36 

Thus, in order to have a non trivial SIC, a description 
in terms of localized orbitals must be adopted, which 
causes violation of the translational invariance and con- 
sequently the inapplicability of the Bloch theorem. Fur- 
thermore, the SIC-LSDA eigenstates are not invariant 
under unitary rotations within the subspace of the oc- 
cupied orbitals, and the SIC-LSDA solutions strictly de- 
pend on the assumption of choosing each orbital as ex- 
tended (thus self-interaction free) or localized (therefore 
subject to SIC). We will show an example of this in Sec- 
tion IV D. 

To our knowledge the implementation of a fully self- 
consistent SIC-LSDA approach for extended systems 
was pioneered by Svane and co-workers. 18 They car- 
ried out a series of applications to a remarkable range 
of materials including the family of transition-metal 
monoxides, 18,2 the high-T c superconductor parent com- 
pounds La 2 Cu04 19 ' 20 and YBa 2 Cu 3 06, 27 the rare-earth 



materials 7-Ce and a-Ce 22 ' 23 ' 25 , and Yb 26 and Pu 28 
monopnictides and monochalcogcnides, obtaining sys- 
tematic improvements over the LSDA results. Other 
SIC-LSDA applications to transition-metal oxides, 21,33 
implemented within different computational methodolo- 
gies, are also present in the literature. 

The major drawbacks of these SIC-LSDA implemen- 
tations are the rather complicated formulation with re- 
spect to the LSDA and the increase of computing cost 
which makes the SIC-LSDA almost impractical for large 
systems. One reason for this increase is the use of big su- 
percells needed to describe the localized orbitals (e.g. ~ 
500 atoms for bulk MnO 19 ' 33 ). Indeed, in the works pre- 
viously cited, the SIC-LSDA is implemented using a basis 
of linear muffin tin orbitals within the atomic sphere ap- 
proximation (LMTO-ASA), whereas, to our knowledge, 
there are no examples of implementations using the more 
expensive plane-wave basis and pseudopotentials. 

An important step towards a practical (albeit approx- 
imate) expression for the SIC-LDA was accomplished by 
Vogel and co-workers 29 by incorporating the atomic SIC 
within the non-local pseudopotential projectors (SIC- 
PP) generated from the free atom. The idea under- 
lying this approach is that the SI potential of a local- 
ized electron in the crystal can be well approximated by 
the SI which the same electron experiences in the free 
atom. The SIC-PP turned out to be quite efficient for 
describing the properties of some highly ionic compounds 
with atomic-like, poorly hybridized bands, such as II- VI 
semiconductors, 29 III-V nitrides, 30 and silver halides. 31 
For these materials, the energy band structures calcu- 
lated with the SIC-PP show a much better agreement 
with photoemission spectroscopy measurements than the 
LDA band structure calculations. 

Our pscudo-SIC can be considered a generalization of 
the SIC-PP approach. It is still based on the idea of 
replacing the SIC potential with a non-local projector, 
but now the projector depends on the orbital occupation 
numbers calculated self-consistently within the crystal. 
The details are described in the next section. 



III. FORMULATION OF THE PSEUDO-SIC 
A. Kohn-Sham equations within pseudo-SIC 

Although the calculations reported in this paper have 
been performed within the USPP method, 34 here we de- 
scribe the pseudo-SIC formalism adapted to NCPP. This 
allows a simpler formulation and an easier understanding 
of the logic behind this approach. (The generalization to 
USPP is given in the Appendix.) 

Within pseudo-SIC, the SIC potential is cast in terms 
of a non-local projector, which resembles the non-local 
part of the pseudopotential: 

VSxc[n,m] - V a HXC \nM - £ M V hxcK] (3H (2) 
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Here n and m are the periodic charge and magneti- 
zation densities of the crystal and the SI is written in 
terms of atomic quantities only: i = nn), R,] is a cu- 
mulative index for angular momentum quantum numbers 
and atomic coordinates, 3^ are projector functions (e.g. 
spherical harmonics), and nf the charge densities of the 
(pseudo) atomic orbitals fc: 



<(r) 



(3) 



where pf are orbital occupation numbers. Through 
Eq.2, the Bloch wavefunctions are projected onto the ba- 
sis of the atomic orbital charges nf . For each projection, 
the Bloch state is corrected by an amount corresponding 
to the atomic SIC potential Vjj XC [nf]. Thus the SIC is 
applied without really introducing a dependence of the 
KS Hamiltonian on the individual Bloch wavefunctions, 
and the difficulties of the SIC-LSDA approach are over- 
come. 

The presence of pf in Eq. 3 is a major novelty of our 
approach. In Rcf. 29 the occupation numbers are implic- 
itly set to 1, i.e. the SIC potentials are generated from 
fully occupied atomic states. Clearly, this cannot be a 
correct choice in general, since the occupation numbers 
may become fractions whenever hybridization, degener- 
acy, or spin-polarization effects arise. Furthermore, mov- 
ing from the free atom to the crystal, the pf can change 
a great deal. For instance atomic orbitals which are fully 
occupied in the atomic ground state may become con- 
duction states in the crystal (e.g. Zn 4s in ZnO). 

Thus the pf must be allowed to be fractional and 
must be recalculated self-consistently within the crystal. 
Within a plane-wave basis set it is straightforward to cal- 
culate pf as atomic orbital projections onto the manifold 
of the occupied Bloch states: 



p1 = £ f& «kl&> <&Kk>- 



(4) 
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These quantities are analogous to the local orbital oc- 
cupations calculated within the FLAPW implementation 
of the LDA+U. 52 ' 53 In the limit of an isolated atom 
the pf recover the atomic values, whereas in the case 
of hybridized bonds, they rescale the atomic SIC by the 
amount of charge that actually occupies the atomic or- 
bital. This fact has a fundamental consequence: since in 
most cases empty bands have dominant characters from 
orbitals whose occupation numbers are close to zero, the 
SIC potential will not affect these bands. This is consis- 
tent with the general assumption that conduction states 
are itinerant, thus Sl-free. (In fact, even if the conduc- 
tion states are not itinerant, it is theoretically justified to 
apply the SIC only to the occupied states, since these are 
the states which actually see their own charge. 48 ) Notice 
that partially occupied bands can be treated on the same 
footing as filled bands, thus the scheme can be applied 
to both insulators and metals. 

Following the suggestion of Ref. 29, for the purpose of 
numerical efficiency we recast the pseudo-SIC as a fully 
non-local, Klcinman-Bylander projector: 



sic 
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where 



and 



7f« = VZ xc [nf(r)}Mr) 



(5) 



(6) 
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Cf = (<t>i\Vg xc [nf}\cf> i ). 
The pseudo-SIC KS equations are: 

-V 2 + V PP + VSxc - V&cl ICk) = <k Kk>, (8) 



where Vpp is the pseudopotential projector, and ef k 
are the KS eigenvalues. 

The recalculation of V^ xc [nf] at each iteration of the 
self-consistency for each atom and angular component 
would result in a major increase of computing cost. A 
large saving of time can be achieved by assuming a linear 
dependence of the SI potential on the occupation num- 
bers: 



VSxcK] = P°VZ xc [nf;pf = l] 



(9) 



so that Vff XC [nf; 1] (i.e. the SI potential for the fully- 
occupied orbital) is set in the inizialization, and only 
the pf needs to be updated during the self-consistency. 
Eq.9 is exact for the Hartree term, which is the domi- 
nant contribution for large occupation numbers, whereas 

1 /3 

it introduces a non-linearity error of0(pt /0 -pi) in the 
exchange-correlation part. 17 

The time saving provided by the linear scaling ar- 
gument is instrumental for calculating structural relax- 
ations and electronic properties of large-sized systems. 
Furthermore, the assumption of linear scaling allows us 
to introduce relaxation effects in a very simple way into 
the pseudo-SIC scheme. Indeed if an electron state is lo- 
calized its energy will change with the orbital occupation. 
Thus, in order to compare the calculated eigenvalue with 
the photoemission spectroscopy data (i.e. with the elec- 
tron removal energy) the effects of the electron relaxation 
must be subtracted out of the one-electron potential. In 
DFT the energy required to remove a fraction p of an 
electron from a one-electron localized state is: 37-39 



AE(p) = E{p) - E(0) 



dte(t), 



(10) 



t=o 



where e is the corresponding KS eigenvalue. The lead- 
ing dependence of the LSDA KS potential (and eigenval- 
ues) on the orbital occupations is indeed contained in its 
SI part. Thus, if 5e is the SI part of the LDA eigenvalue, 
within linear scaling we have AE(p) = 1/2 p 2 Se(l). It 
follows that the SIC relaxation energy, 1/2 p 2 5e(l), can 
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be subtracted out by rescaling the SIC potential as fol- 
lows: 



V, 



HXCl 



\v HXC K 



(11) 



Through Eq.ll we directly incorporate in the pseudo- 
SIC KS Equations the electron removal energy due to 
the SI contribution. Wc point out that the eigenvalue 
relaxation is very important in order to match the pho- 
toemission spectroscopy results. Discarding this contri- 
bution, the SIC eigenvalues would strongly overestimate 
the electron removal energies. 

Notice that, if the KS eigenvalues did not depend on 
the occupation numbers, according to Eq.10 they would 
be equal to the electron removal energies. This is in fact 
the case in Hartree-Fock calculations and in any other 
theory which obeys Koopman's theorem. This property 
does not hold in LDA or any LDA-related scheme (such 
as GGA, SIC, etc.) where the potential explicitly de- 
pends on the orbital occupation numbers. 

As is customary in atomic calculations, we assume the 
radial approximation for the atomic orbital charges, so 
that the SIC projectors can be written: 

li, mi A r ) = \pl mi , v V^ xc [nl v {r)-,l]ci>l muV {v), (12) 

where v labels the atom type, and nf„(r) is the ra- 
dial pseudo-charge density of orbital (I, mi) (in radial ap- 
proximation this does not depend on m;). Finally, the 
normalization coefficients are: 
1 „ 



In summary, except for the electron occupation num- 
bers, all the other atomic-like ingredients can be im- 
ported from the code used for the pseudopotential gen- 
eration, and the Kleinman-Bylander projectors for the 
SIC are set during the initialization process. Since the 
calculation of the pf m v through Eq.4 is not very de- 
manding, the global computational cost of the pseudo- 
SIC for each self-consistent iteration is roughly equal to 
that of the basic LSDA. However within pseudo-SIC the 
number of self-consistent iterations required to reach the 
self-consistency can be larger than in LSDA due to oscil- 
lations in the values of the occupation numbers. 



B. Total Energy within pseudo-SIC 

In the previous section we estabilished the form of the 
pseudo-SIC KS Equations. Here we formulate a suitable 
expression for the total energy functional. We point out 
that, within our scheme, a physically meaningful energy 
functional which is also related to Eqs. 8 by a variational 
principle is not available. (It is, by construction, within 
LSDA and SIC-LSDA.) 

In SIC-LSDA the energy functional is: 17 



E S ic[n,m} = E[n,m] 



E E hxcK 



(14) 
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where E[n, m] is the LSDA energy functional and 
-Ehxc [nf] the Hartree exchange-correlation energy of the 
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EhxcK] = dvn^v) - V H K(r)] + £xcK(* 



(15) 



where Vh is the Hartree potential and £xc is t ne local exchange-correlation energy density. For the pseudo-SIC 
total energy we adopt the same expression as Eqs. 14 and 15, with the orbital charges nf given by Eqs. 3 and 4. 
Notice that, 



SE H xcK] 
Spf 



= Cf. 



(16) 



Eq.16 represents the Janak thoerem 37 ' 39 applied to the SIC contribution since Cf (see Eq.7) is the SI part of the 
atomic eigenvalue. In terms of the pseudo-SIC eigenvalues (Eq. 8) the total energy can be rewritten: 

Esicln, m] = Y, <k ~ E / dr n ' » V HXC [n(v),m(r)} + E HXC [n, m] + E ion 



+ \Y. fnU «k I VSIC I <k) - E E **XC [<] , (17) 



2 , 

nk., a 



•5 



where Ei on is the usual Ewald term. Y_\ a Enxcinf] 
produces a gentle modification to the LSDA energy func- 
tional whereas the fifth term in Eq.17 is a strongly vary- 
ing contribution which compensates the same contribu- 
tion present in the pseudo-SIC eigenvalues. Without this 
compensation, Eq.17 would give very inaccurate total en- 
ergies which would be unphysically far from the LSDA 
values. A numerical example of this behavior will be 
shown in Section IV C. 



IV. RESULTS 

A. Technical details 

In this work we compare results from our plane 
wave USPP 34 pseudo-SIC implementation with results 
from conventional LSDA USPP method. The local 
exchange-correlation energy functional is modeled using 
the the Perdew-Zunger interpolation formula. 17 The use 
of USPP 34 allows us to obtain well -converged results with 
moderate cut-off energies (35 Ry for MnO and YMn03, 
40 Ry for ZnO and 45 Ry for GaN). In order to have 
highly transferable USPP, two projectors per angular 
channel are included for all the atoms. The 3d 10 electrons 
in ZnO and GaN are treated as valence states, while the 
semicore Y s and p electrons are placed in the valence for 
the LSDA calculations and in the core for the pseudo-SIC 
calculations (the reason will be explained in the follow- 
ing section). For total energy calculations we use up to 
8x8x8 grids of special k-points. 40 



B. Atomic ingredients for the pseudo-SIC 

The atomic quantities necessary to build the pseudo- 
SIC hamiltonian are the pseudo atomic orbitals 4>i^ mu v, 
required in the Kleinman-Bylander projector and for the 
calculation of pf , and the SIC potentials V£ xc [nf v (r) ; 1] 
for the respective pseudo orbital charges calculated at full 
electron occupation (we use pseudopotential and not all- 
electron atomic functions for obvious reasons of smooth- 
ness.) 

To illustrate the impact of the atomic SIC, in Figure 
1 we compare the Zn s, p, and d (unpolarizcd) pseu- 
dopotcntials (V;) to the same quantities minus the cor- 
responding corrections AV^ 10 = VHxc[ n iyA]- As ex- 
pected, the SIC makes the electron potentials more at- 
tractive and recovers the physical long-range limit -2/r 
(inRy). 

Note that, before transfering VuxcV^iy, 1] to the crys- 
tal, the long-range tails must be cut off, since the SIC 
should act as a local correction: each SIC potential must 
be directly applied only to the electron states localized 
on the same atom, otherwise the overlapping Coulomb 
tails would give rise to an unphysical SIC overestimation. 
In Figure 2 we report Vuxc[niy 1] (times r) for both 



pseudo and all-electron orbital charges (the latter show 
the cusps corresponding to the wavefunction nodes), and 
the pseudo orbital charges p\. Since only the products 
Vhxc[ u Im] Pi contribute in the Kleinmann-Bylander con- 
struction, the SIC potentials can be cut off as soon as the 
orbital charges vanish, without losing accuracy. 

2 r , 
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FIG. 1. Zn pseudopotentials (V;) and the corresponding 
atomic SIC (AV, S/C ) calculated at full orbital occupation. 
The correction lowers the one-electron potential and recovers 
the physical long-range limit -2/r which is missing in LDA. 

It is important to notice that Vff XC [nf v ; 1] is not very 
sensitive to the atomic valence configuration from which 
it is actually generated. For example the integrated value 

°f ^Hxc[ n 3 d ' 1] ' <^3d (i- c - tne a t om i c SI) changes by only 
- 1% if calculated in the 3d 10 4s 2 , the 3d 10 4s 1 , or the 
3d 10 4s° configuration. In other words, despite being va- 
lence dependent, V^ xc [niy, 1] is fairly transferable. 

Notice also that the pseudo-SIC does not interfere with 
the pseudopotential construction, and the same pseu- 
dopotentials used for the LSDA calculations can be used 
for the pseudo-SIC as well, due to the presence of the 
occupation numbers. Indeed, in the limit pi—0 the SIC 
vanishes, and the hamiltonian must continuously recover 
its LSDA value. 

However for pseudo-SIC calculations we prefer to build 
pseudopotentials which take into account the SIC for the 
core states. Although this correction does not signifi- 
cantly modify the valence properties, 41 there are reasons 
related to pseudopotential transferability which suggest 
that the SIC should be included in the core states. First, 
the SIC could push strongly localized valence states (e.g. 
the semicore states) down in energy and unphysically 
close to the core states if the latter are not shifted down 
consistently. Second, the application of SIC to the semi- 
core states may better justify their inclusion in the core 
even if in LSDA they need to be treated as valence states. 
This is the case for the Y 4s and 4p electrons, and the 
Ga 3d electrons. 
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FIG. 2. Zn atomic SIC for psoudo (solid line) and 
all-electron (dashed line) orbital charges. The pseudo orbital 
charges (dotted lines) are plotted on the positive y-axis. The 
pseudo SIC are cut-off at a diagnostic radius where the orbital 
charges vanish. 

A third reason is even more fundamental: increasing 
the binding energy and the space localization of the core 
states makes them less sensitive to charge relaxation or 
chemical activity. This reinforces the hypothesis that 
core states are insensitive to the chemical environment, 
which is at the basis of the pseudopotential approxima- 
tion. 

The transferability of the atomic SIC to the extended 
system can be tested with the same procedure used for 
the pseudopotentials: the atomic pscudo-SIC hamilto- 
nian should recover the all-electron, fully self-consistent 
SIC-LSDA eigenvalues. We generally find that pseudo- 
SIC and SIC-LSDA atomic eigenvalues agree within ~ 
1% 

Finally we stress a fundamental point: our proce- 
dure is not in any way restricted to the use of pseu- 
dopotentials or any other particular choice of basis func- 
tions. For example it can be equally well implemented 
within the context of all-electron methods. The ad- 
vantage of using the pseudopotential approach relies on 
the fact that the SIC projectors can be easily assembed 
from quantities (such as the projections of Bloch states 
onto atomic orbitals) which are usually calculated in any 
pseudopotcntial-based computing code. 



C. Wide-gap semiconductors 

We begin our series of pseudo-SIC applications with 
two of the most prototypical wide-gap semiconductors, 
wurtzite ZnO and GaN. These materials are ideally 
suited as test cases since the LDA predictions for the 



band energies show some easily recognizable discrepan- 
cies with the abundantly available photocmission spec- 
troscopy data. Furthermore, these compounds are of 
technological interest due to their applications in opto- 
electronic and piezoelectric devices. Nowadays the elec- 
tric polarization and related properties can be efficiently 
calculated by first-principles techniques. 42 Thus a theory 
able to repair the LDA description of the electronic struc- 
ture will also be important for the accurate evaluation of 
dielectric and piezoelectric response. 

Since ZnO and GaN show similar structural and elec- 
tronic behavior, we describe them together. In Table I we 
report our calculations within LDA and pseudo-SIC for 
the equilibrium structural parameters. These have been 
evaluated by energy minimization within the full space of 
parameters a, c and u (u is the anion-cation distance in 
units of c). In general, the SIC volume is larger than the 
LDA volume, since the stronger electron localization due 
to the SIC enhances the electron screening and reduces 
the electron-ion interaction and the charge hybridization. 
This is generally a favorable correction since it is well 
known that within LDA (or LSDA) the lattice parame- 
ters are underestimated compared with the experimental 
values. 

TABLE I. Equilibrium lattice parameters (a,c,u), band gap 
(E s ) and average d-band energy (E d ) for ZnO and GaN. All 
the distances are in Bohr except the anion-cation distance, u, 
which is in units of c. All the energies are in eV. For ZnO the 
experimental values are from Ref. 43. 





LDA 


pseudo-SIC 


Expt. 


ZnO 


a 


6.12 


6.17 


6.16 


c 


9.88 


9.79 


9.84 


u 


0.378 


0.384 


0.382 


E 9 


0.94 


3.70 


3.4 


E d 


-5.3 


-7.5 


-7.8 


GaN 


a 


6.03 


6.045 


6.03 4 ' 


c 


9.80 


9.81 


9.80 47 


u 


0.377 


0.378 


0.375 47 


E 9 


2.16 


4.26 


3.5 43 


E d 


-13.8 


-18.1 


-17.1 49 
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In the case of ZnO and GaN, however, our LDA val- 
ues are already in excellent agreement (within less than 
~ 1%) with the experiments and no correction would be 
required. Happily, the pseudo-SIC calculations do not 
overcorrect the LDA values, (the reason will be explained 
later on the basis of the band structure results). This is 
an important aspect since typical beyond-LDA method- 
ologies, also not derived by variational principles, often 
improve the LDA description of the electron excitation 
spectra but are less accurate than LDA in the determi- 
nation of the equilibrium structure. 

ZnO LDA 
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FIG. 3. LDA band structure for wurtzite ZnO. The orbital 
character for each group of bands is indicated. 



In Figure 3 our calculated LDA energy band struc- 
ture of wurtzite ZnO is shown. For such strongly ionic 
compounds each group of bands can be clearly labeled 
according to a single, dominant orbital character as indi- 
cated in the Figure. 



With respect to the photocmission spectroscopy re- 
sults, the 3d bands calculated within LDA are too high 
in energy and overlap with the sp 3 valence band mani- 
fold. This produces a spurious p-d hybridization which 
shrinks the energy range of the sp 3 bands. Furthermore, 
there is the notorious problem of the fundamental energy 
gap which is underestimated by ~ 40 % (see Table 3). 



ZnO pseudo-SIC 




FIG. 4. Pseudo-SIC band structure for wurtzite ZnO 



All these undesirable features are largely overcome in 
our calculated pseudo-SIC band structure (Fig. 4). The 
largest SIC effect is a downshift of the d band energies, 
which are now placed ~ 3 eV below the center of the 
sp 3 -band manifold, in excellent agreement with the ex- 
periments. Accordingly, the sp 3 band manifold is ~ 1 eV 
broader than in LDA. 

Furthermore, the bands with O p orbital character, 
which are the main contributors to the valence band top 
(VBT), are also corrected for an amount of SI corre- 
sponding to their average electron occupation (~ 0.8 for 
the O p orbitals). In contrast the conduction band bot- 
tom (CBB), which is mainly Zn s in character, is almost 
unchanged, since the Zn s orbital occupation is ~ 0.1. As 
a consequence the pseudo-SIC energy gap opens up to a 
value fairly close to the experimental band gap. 

Our results show a substantive agreement with the cal- 
culated SIC ZnO band structure of Ref. 29 (they obtain 
E g = 3.5 eV). However the two methods act in quite dif- 
ferent ways: in the approach of Ref. 29 all the bands 
(occupied and unoccupied) are corrected by the atomic 
SIC potential of the corresponding fully occupied orbital 
charges. Thus in Ref. 29 the increase of the energy gap 
with respect to the LDA can be roughly quantified as 
the difference between the SIC of the fully occupied Zn 
4s and O 2p atomic eigenvalues. Instead, in our scheme 
only the occupied bands are corrected. We belive that 
this is a more conceptually sound approach for the rea- 
sons explained in Section III A. 

In Figure 5 we show the LDA band structure of 
wurtzite GaN. As is usual in LDA calculations, the Ga d 
band energies fall within the same energy range as the N s 
bands, in disagreement with the experiments which place 
the d bands ~ 3 eV below the N s bands. 49 This gives 
rise to a spurious s-d hybridization which causes an in- 
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crease of lattice constant due to the closed-shell repulsion 
enhanced by the resonances between states on neighbor- 



GaN pseudo-SIC 



mg sites 



44-46 



(i.e. N s and Ga d). In our calculation 



this effect compensates the LDA tendency of underesti- 
mating the lattice parameters, and causes a better-than- 
usual agreement with the experimental data (Tab. I). 
However, the description of the band structure appears 
grossly inappropriate when compared to the X-ray pho- 
toelectron spectra. 49 The s-d hybridization splits the N s 
bands into two sections, one placed above and the other 
below the d bands. The s-d manifold is positioned only 
~ 6 eV below the bottom of the valence sp 3 manifold, 
and the energy gap is 40% smaller than the experimental 
value. 



GaN LDA 
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5. Band structure of wurtzite GaN calculated in LDA. 



The pseudo-SIC corrects in large part the LDA de- 
scription. In Figure 6 we show the GaN band structure 
calculated within pseudo-SIC. Here the Ga d bands (E^ 
~ -18 eV) are correctly located ~ 3 eV below the center 
of the N 2s bands, thus the spurious s-d hybridization 
is avoided and the closed-shell repulsion reduced. This 
explains why the equilibrium lattice parameters calcu- 
lated within pseudo-SIC are close to the LDA values al- 
though the pseudo-SIC is generally expected to enhance 
the electron-ion screening. 

Also, the binding energies of Ga d and N s states are 
increased by ~ 4 eV and 2 eV respectively, thus cor- 
recting in large part the faults of the LDA description. 
However within pseudo-SIC the energy gap is somewhat 
overcorrected (it is ~ 0.9 eV larger than the experimental 
value). 




FIG. 
pseudo- 
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Band structure of wurtzite GaN calculated in 



SIC. 
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FIG. 7. Total energy values for zincblende GaN as a func- 
tion of the lattice constant ao. The meaning of the curves is 
explained in the text. 

To conclude this Section, in Figure 7 we analyze the 
behavior of the pseudo-SIC total energy (Eq.17) for 
zincblende GaN as a function of the lattice constant in 
the region around its experimental value ao=8.505 Bohr. 
Ei is the functional given by the first four terms of Eq.17. 
Clearly, its minimum value is far out of the physically 
meaningful region. E2 is obtained by adding to Ei the 
fifth term of Eq.17, thus E 2 equals the LSDA energy 
functional (except for the difference in the wavefunc- 
tions) . This functional has a minimum value in excellent 
agreement with the experiments. Finally, Eg/c is E2 + 
J2i a ^hxc [nf] ■ This last term is almost independent of 
the lattice parameter, since it only changes through the 







pi, which are rather local quantities. As a consequence, 
Egjc and E2 have similar behavior. 



MnO LSDA 



D. Transition-metal oxides. 



The LSDA misrepresentations of the electronic proper- 
ties of the transition-metal oxides originate mainly from 
the description of the d electron states. At variance with 
the non-magnetic semiconductors considered in the pre- 
vious section, in the transition-metal oxides the d states 
lie higher in energy and closer to the fundamental gap, of- 
ten overlapping with the oxygen p states. Since in LSDA 
the d electron binding energies are severely underesti- 
mated, the d character is generally dominant at the VBT, 
so these compounds are described as Mott-Hubbard in- 
sulators. This is in striking contrast with the photoemis- 
sion spectroscopy data which describe these materials as 
charge-transfer insulators (or eventually in the interme- 
diate charge-transfer Mott-Hubbard regime) with a ma- 
jority or even dominant O p character at the VBT. 7 



Furthermore the LSDA energy gaps for these materi- 
als are even more severely underestimated than those of 
the III-V or II- VI semiconductors and in some cases the 
gap can be vanishing. This is because the major driving 
force leading to the formation of a gap separating bands 
of equal orbital character should be the on-site Coulomb 
energy U, but in LSDA the energy gap can only open due 
to Hund's rule and the crystal field splitting. These are 
both of order 1 eV, i.e. almost one order of magnitude 
smaller than U, and therefore can be easily overcome by 
the the band dispersion. 

Finally, since the LSDA tends to broaden the space 
distribution of the electron charge and emphasize the in- 
teratomic charge hybridizations, the local magnetic mo- 
ments are generally underestimated with respect to the 
experimental values. 
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FIG. 8. LSDA band structure of MnO. The orbital charac- 
ter is indicated for each band. The orbitals are expressed in 
hexagonal coordinates, i.e. z is parallel to the [111] direction. 



As applications of the pseudo-SIC to magnetic com- 
pounds we choose two of the most widely studied 
transition-metal monoxides, MnO and NiO. For these 
materials we can compare our results with both ex- 
perimental data and the output of some common 
methodologies used for repairing the faults of the 
LSDA description. 5 ' 7 These schemes include the SIC- 
LSDA, 18 ' 24,33 the LDA+U, 32 ' 24 ' 51-53 and the model 

GW 54-57 

The ground state of MnO and NiO is characterized 
by A-type antiferromagnetic ordering, which consists of 
(111) ferromagnetic layers of Mn (or Ni) with alternating 
spin directions, while the oxygens are non-magnetic. The 
symmetry is rhombohcdral, with 4 atoms per primitive 
unit cell. In our calculations we fix the lattice constant 
to its experimental value (ao = 8.37 Bohr and 7.93 Bohr 
for MnO and NiO, respectively). 

In Figure 8 the LSDA band structure for MnO is 
shown. Notice that the orbital characters are referred 
not to the cubic, but to the hexagonal coordinates, i.e. x 
and y lie on the hexagonal (111) plane, and z is parallel 
to the [111] axis. Also, within the hexagonal (or rhom- 
bohcdral) crystal field, the 3d orbitals are split into two 
doublets (d xv , d x 2_ y 2) and (d xz , d yz ), and one singlet 

The nominal ionic configuration Mn 2+ 2_ suggests a 
Hund's rule-induced splitting between filled spin up and 
empty spin down 3d bands, resulting in an energy gap. 
From the analysis of the band character we see that the 
VBT is composed of 60% (d^, d[ 2 _ y2 ) and (d^, djj 
orbitals, with the remaining 40% coming evenly from the 
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(j>x, Py) doublets of the two equivalent oxygens. In con- 
trast, the conduction band bottom (CBB) consists of ~ 
80% d^ 2 centered on the same Mn, which is equivalent 
to the d^ 2 orbital localized on the other Mn. Thus, from 
the LSDA calculation, a picture of the MnO as a small- 
gap Mott-Hubbard insulator emerges, which is not what 
should be expected according to the experiments. 

Furthermore, while the value of the Mn magnetic mo- 
ments is satisfactory, the size of the LSDA energy gap 
severely underestimates the measured value (see Table 
II). Also, notice that the occupied 3d bands are rather 
fiat and well separated by a ~ 1 eV gap from the under- 
lying O 2p band manifold which is ~ 6 eV wide. 

In Figure 9 we see how the application of the pseudo- 
SIC modifies the results. First of all, there is no longer a 
gap between 3d and 2p bands. The SIC causes a down- 
ward shift of the Mn 3d bands with respect to the more 
dispersed O 2p bands, thus the 3d character is now spread 
across the whole (8 eV wide) valence band manifold, and 
is heavily mixed with the O 2p character. 

This shift significantly increases the energy gap which 
is now well within the experimental uncertainty (see Ta- 
ble II). It also changes the VBT character, which is now 
composed of 40% 3d doublets and 60% 2p orbitals. 85% 
of the CBB still comes from the d^ 2 orbital. 

Thus the agreement with experiments, which locate 
the MnO in the intermediate charge-transfer Mott- 
Hubbard regime, is restored. Furthermore, the enhanced 
on-site localization of the 3d orbitals due to the SIC in- 
creases the magnetic moment which is now in excellent 
agreement with the experimental value. 

TABLE II. Magnetic moments, M, and energy band gaps, 
E 9 , of MnO and NiO. The upper part shows our results calcu- 
lated within LSDA and pseudo-SIC, in comparison with the 
experimental values. The lower part shows results of other be- 
yond-LSDA calculations (results in parentheses are explained 
in the text). 

MnO NiO 





E 9 (eV) 


M(/is) 


E 9 (eV) 




LSDA 


0.92 


4.42 


0.4 


1.11 


pseudo-SIC 


3.98 


4,71 


3.89 


1.77 


Expt. 


3.8-4.2 60 


4.79, 58 4.58 s9 


4.0, B1 4.3 62 


1.77 s8 1.90 s9 


SIC-LSDA 18 ' 21 


3.98 


4.49 


2.54 


1.53 


SIC-LSDA 33 


6.5(3.4) 


4.7(4.7) 


5.6(2.8) 


1.7(1.5) 


LDA+U 32 


3.5 


4.61 


3.1 


1.59 


LDA+U 53 






4.1(2.8) 


1.83(1.73) 


LDA+U 52 






3.38 


1.69 


GW 55 


4.2 


4.52 


3.7 


1.56 
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MnO pseudo-SIC 




FIG. 9. Pseudo-SIC band structure of MnO. The orbital 
character is indicated for each band. 



The LSDA description of NiO is even worse than that 
of MnO. Indeed, the Ni 2+ ion has nominally eight elec- 
trons in the occupied bands, thus the energy gap within 
the d band manifold can only be opened by the weak 
crystal field splitting. 

In Figure 10 we show the band structure of NiO ob- 
tained within LSDA. The fundamental gap occurs be- 
tween bands derived from the singlet d 2 and the doublet 
(d^ z ,d^ z ). In particular, the VBT is almost purely (90 

%) d^ 2 in character, and no hybridization with oxygens is 
present. The gap value (see Table II) is only 0.4 eV, and 
in fact some LSDA calculations even describe NiO as a 
metal. 33 The majority d bands are rather flat and, as in 
MnO, separated by a ~ 1 eV wide gap from the p bands. 
Finally the magnetic moments are underestimated by ~ 
40%. Thus, according to the LSDA calculations, NiO is 
a small-gap Mott-Hubbard insulator, whereas the exper- 
iments describe NiO as a wide-gap charge-transfer insu- 
lator. 

The inclusion of SIC (Figure 11) greatly improves the 
NiO description. As in MnO, the eight occupied d bands 
are shifted down in energy and strongly hybridized with 
the O p bands. The VBT becomes a mixture of d and 
p character, with a predominance of the latter. For ex- 
ample at r the VBT singlet is almost purely p 2 , whereas 
at K the twofold degenerate VBT comes from the hy- 
bridization between (p x , p y ), which contribute 80%, and 
(d], z , dj z ), which contribute 20%. Thus, it is a pd-zr- 
hybridized band. The same pd7r hybridization character- 
izes the CBB doublet, but the percentages of the compo- 
sition are reversed: 90% of the charge comes from (d* z> 
d* z ), and only 10% from (p x , p y ). Thus, according to the 
pseudo-SIC, NiO is a charge-transfer insulator. Further- 



more, the calculated energy gap and magnetic moments 
are in good agreement with experiments (see Table II). 
As is evident from the band structure, the Ni magnetic 
moment, M=1.77 /is, comes from the spin-polarization 
of the degenerate orbitals d xz and d yz , each of them car- 
rying M/2 magnetization. Finally, the bands at the bot- 
tom of the occupied p-d manifold are 90% from the dou- 
blet (dl y , d^ 2 _ 2/2 )- These orbitals lie on the hexagonal 
plane and are entirely localized since they do not point 
towards the oxygens. As a result their corresponding 
bands experience the largest SIC. 



NiO LSDA 
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FIG. 10. LSDA band structure of NiO. The orbital charac- 
ter is indicated for each band. 

A number of beyond-LSDA approaches have been at- 
tempted in the past and tried on the series of transition- 
metal oxides. In Table II we listed the results of some of 
these calculations. 

Whithin full SIC-LSDA, the results for the transition- 
metal oxides are critically sensitive to the choice of or- 
bitals taken as localized. In Ref. 33 (a LMTO-ASA cal- 
culation) the SIC-LSDA band structures are calculated 
in two ways: The first calculation assumes that 3d and 
the 2p orbitals are both localized and therefore affected 
by the SI, and the second is performed with only the 3d 
orbitals taken as localized. (In the Table we report the 
results of both the calculations; the d-only SIC-LSDA re- 
sults are in parentheses). It is shown that the first option 
gives a band structure in better agreement with photoe- 
mission spectroscopy, with the exception of the band gap 
which is strongly overestimated. Instead, the choice of 
correcting only the d orbitals gives a better energy gap 
but also leads to a too large downshift of the d bands 
with respect to the p bands (e.g. for MnO and NiO the 
d bands lie ~ 5-6 eV below the center of the p bands 33 ). 
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As a consequence, the calculation exaggerates the charge- 
transfer character of the band gap. 



NiO pseudo-SIC 




FIG. 11. Pseudo-SIC band structure of NiO. The orbital 
character is indicated for each band. 



The strategy of correcting both d and p orbitals repre- 
sents our preferred point of view, 50 and is consistent with 
the pseudo-SIC which assumes that all the electron states 
should be corrected, without any a-priori discrimination. 
Indeed, the pseudo-SIC band structures are similar to 
those obtained in Ref. 33 with the first choice of orbitals, 
with the important exception that in our pseudo-SIC 
the energy gap is in much better agreement with experi- 
ments. We speculate that the reason for this difference is 
the absence in SIC-LSDA of the relaxation energy con- 
tribution, which shifts the occupied d and p bands up in 
energy. 

The LDA+U corrects only the d states, and it shows 
similarities with the d-only SIC-LSDA in describing the 
properties of the transition-metal oxides. 32,24,51,53 Indeed 
in both these approaches, the transition-metal oxides are 
all described as charge-transfer insulators. However in 
LDA+U the amplitude of the d band shift is critically 
controlled by the parameter U. For example, in Ref. 53 
(a projector augmented-wave (PAW) calculation) it is 
shown that for NiO the widely used value U=8 eV gives 
a good energy gap but a poor description of the occupied 
band manifold due to the exaggerated downward shift 
of the d bands. In contrast, the value U=5 eV gives a 
smaller energy gap, but optical properties and magnetic 
moments in better agreement with experiments. (Results 
for U=5 eV are reported in parentheses in Table II). 

Finally, we report the results of the model GW, 54 
an approach radically different from both the SIC or 
the LDA+U schemes, based on a model self-energy 



correction to the LSD A KS equations. The model 
GW is capable of giving an accurate description of 
both transition-metal oxides 55,56 and other non-magnetic 
semiconductors, 54,55,57 thus it should probably be con- 
sidered as the most reliable reference. The fact that the 
pseudo-SIC values for the energy gaps are close to the re- 
sults obtained within GW is an indication that, at least 
for this family of compounds, the SI, and not many-body 
effects, really represents the main source of the LSDA er- 
ror. 



E. Hexagonal manganites 



As a final application of the pseudo-SIC, we con- 
sider hexagonal YMn03, which has recently attracted 
the attention of both the experimental 63-66 and the 
theoretical 67,11 communities since it shows the uncom- 
mon characteristic of being both magnetically and fer- 
roelectrically ordered within the same bulk phase. The 
study of YMn03 is motivated by the possible applica- 
tions of ferroelectromagnetic materials as building blocks 
for spintronic devices 68 and, from a more fundamental 
standpoint, by the necessity to understand the interac- 
tion between magnetic and ferroelectric polarization and 
the conditions which favor this coexistence. 69,11 

YMn03 can be grown in both orthorhombic 70 and 
hexagonal 71-73 structures, although the latter is the most 
stable. The orthorhombic phase, typical of manganese 
perovskites (e.g. LaMnOs), is antiferromagnetic but not 
ferroelectric. The hexagonal phase shows a spontaneous 
ferroelectric polarization P ~ 5.5 /j,C/cm 2 parallel to the 
c axis 74,72 below a critical temperature T c =900 K, and 
is A-typc antiferromagnetic 75,76 at temperatures below 
Tjv=80 K. In this section we will consider only hexago- 
nal YMn0 3 . 
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FIG. 12. Structure of YMnC>3 in the paraelectric P63/mmc 
phase. The arrows are placed on the Mn ions and indicate the 
spin-polarization direction (thus the ordering is A-type anti- 
ferromagnetic). Each Mn is surrounded by a bi-pyramidal 
cage of five corner-sharing oxygens. Atoms not connected by 
bonds are Y. 
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FIG. 13. LSDA band structure for the A-type AFM para- 
electric YMn03. The orbital character is given in terms of 
the hexagonal cartesian coordinates. According to this calcu- 
lation the system is a metal. 



The structure of the paraelectric YMnC>3 (P63/mmc 
symmetry) is shown in Figure 12. The Mn ions, sited 
on close-packed hexagonal positions, are surrounded by 
corner-sharing bipyramidal cages of oxygens. The unit 
cell is made of eight lxl hexagonal planes enclosing a 
total of 10 atoms. In the ferroelectric phase (P63CUL sym- 
metry) the MnOs bipyramids are tilted around the axis 
passing through the Mn and parallel to one of the trian- 
gular base sides, thus the hexagonal planes are y/3 x V3 
and the unit cell has 30 atoms. However, for the pur- 
poses of illustrating the effects of the pseudo-SIC, the 
small oxygen rotations of the ferroelectric structure are 
not significant, thus in the following we only consider the 
paraelectric phase. 



The technical features of our calculations are the same 
as those used in Ref. 11, where a detailed study of den- 
sity of states, band structure and orbital charges within 
LSDA can be found. Here we focus on the remarkable 
changes in the band structure of this material produced 
by the pseudo-SIC. 
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FIG. 14. Pseudo-SIC band structure for A-type AFM, 
paraelectric YMnOa. The pseudo-SIC opens a gap between 
empty d bands and filled pd er-hybridized bands. 

An obvious requirement for a ferroelectric material is 
that it must be insulating, but the LSDA calculations de- 
scribes YMnOa as a metal, as can be seen in Figure 13. 
(The hexagonal crystal field splitting and the cartesian 
coordinates are the same as those described previously 
for MnO and NiO.) In hexagonal symmetry the four d 
electrons of the Mn 3+ ion entirely occupy the two or- 
bital doublets (dl y , d} x2 _ y2 ) and (dl z , d) yz ), leaving the 

d , a orbital, which is the highest in energy, empty. This 
ordering causes a magnetic moment on Mn equal to 3.8 
fiB, slightly lower than its nominal value 4 /ib because 
of Mn d-0 p hybridization. 

However the LSDA crystal field splitting is too small 
to open a gap, and one band crosses Ep in the direction 
T-Q-P which is parallel to the k z —0 plane. This band 
comes from a mix of Mn d\ y , d\ 2 _ 2 orbitals and p£, pj, 
orbitals from the oxygens lying in-plane with the Mn (i.e. 
it is a pdfj hybridization). In contrast, the band energies 
are fiat along [0001], as is typical for strongly layered 
compounds. (Notice that according to our LSDA calcu- 
lations the system is metallic even within the ferroelectric 
phase.) 

The metallicity is a fatal shortcoming since it precludes 
the possibility of accessing any ferroelectric properties, 
such as spontaneous polarization, Born effective charges, 
and piezoelectric tensors, (according to our calculations 
the system is still metallic within the ferroelectric struc- 
ture) 

The pseudo-SIC (Figure 14) repairs the fault: a gap, 
E g = 1.40 eV, opens (in excellent agreement with the ex- 
perimental value Eg— 1.47 eV) between the empty d\ 2 
band and the fully occupied pdc band. Thus the pseudo- 
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SIC calculations describe the paraelectric YMnOa as an 
intermediate Mott-Hubbard charge-transfer insulator, in 
agreement with photoemission experiments 73 and previ- 
ous LDA+U calculations. 67 Also, the two filled doublets 
(dl y , d} x2 _ y2 ) and (dj. z , d yz ) are shifted down by ~ 3-4 eV 
with respect to the LSD A values, and are located below 
the O p manifold. 

On the basis of the pseudo-SIC calculation which cor- 
rectly describe the YMnC>3 as an insulator we are now 
able to study the ferroelectric structural displacements 
and other properties related to the electric polarization. 
This will be the subject of future work. 



V. DISCUSSION AND SUMMARY 

In this paper we have proposed and implemented an 
innovative first-principles approach which is able to im- 
prove the LDA and LSDA description of the electronic 
properties for a vast range of compounds, including 
strongly-correlated and magnetic systems, and requires 
only a minor increase of computing cost. Our strategy, 
based on the approximation of the true SIC potential 
with a pseudopotential-like projector, was built primar- 
ily on the work of Rcf. 29, where it was shown that the 
inclusion of the atomic SIC within the pseudopotential 
construction can consistently improve the agreement of 
band structures with photoemission spectra for a range of 
strongly ionic II- VI and III-V insulators. Working on this 
original idea, we defined a pseudo-SIC functional explic- 
itly dependent on the orbital occupation numbers which 
are self-consistently calculated from the Bloch wavefunc- 
tions and represent the natural extension to periodic sys- 
tems of the atomic occupation numbers. They can be 
fractional due to charge hybridization, eigenvalue degen- 
eracy, or Fermi-Dirac distribution. 

We tested the pseudo-SIC on three different classes of 
compounds: non-magnetic semiconductors, magnetic in- 
sulators and ferroelectrics, and found a generally good 
(sometimes very good) agreement with experimental 
data. In particular, the size and the orbital characteri- 
zation of the fundamental energy gap and the local mag- 
netic moments are much better described than in LSDA. 

We expect the same kind of improvement in any system 
where the electron localization plays a major role. Exam- 
ples of systems which could benefit from the pseudo-SIC 
application are innumerable, including bulk systems with 
d and/or f electrons, Jahn- Teller and orbitally ordered 
systems, and non-bulk systems like defects and impu- 
rities, surface and interface states, bulk resonances and 
core states. The SIC may affect a vast range of phenom- 
ena, encompassing, among others, surface reconstruc- 
tions, adsorbtion and diffusion of atoms and molecules 
on surfaces, doping in semiconductors, alloys, and homo- 
and hetero-j unctions. 

Of course, the discrepancies between theoretical and 
experimental results which are not attributable to the SI 



should not be expected to disappear within pseudo-SIC. 
A typical example is the fundamental energy gap of bulk 
Si, whose LDA value is ~ 0.5 eV lower than the experi- 
mental gap, due to genuine self-energy effects which are 
outside the realm of the DFT itself. As a consistency 
test, we calculated the band structure of bulk Si within 
pseudo-SIC. Happily, it comes out to be very similar to 
the LDA band structure. This is a consequence of the 
fact that each band (either occupied or empty) shares 
the same orbital character (sp 3 -hybridization) , thus the 
pseudo-SIC mostly produces a trivial shift of the whole 
band manifold. 

Comparing it to other corrective methodologies avail- 
able in the literature, the pseudo-SIC is closest to the 
spirit of the LDA+U: in both cases the LSDA energy 
functional is augmented by a term depending on the or- 
bital occupation numbers. In the pseudo-SIC the modi- 
fication is milder, and consists of subtracting out a con- 
tribution (i.e. the SI) which is already present in the 
LSDA energy functional. The LDA+U is a more radi- 
cal departure from the LSDA, since it rewrites the whole 
Coulomb electron-electron interaction in terms of local 
orbitals according to the multiband Hubbard expression 
(which is, in itself, self- interaction free). Based on the 
undeniably limited set of results presented in this paper, 
the pseudo-SIC seems to be at least as accurate as the 
LDA+U, but a fair comparison will need a much larger 
body of work. 

Finally, we point out that our method suffers a 
drawback, that is the non-variationality of the energy 
functional. This is especially annoying when force or 
stress calculation is required, since the familiar Hcllman- 
Feynman formulations cannot be applied and additional 
contributions arise due to the first-order change of the 
wavefunctions. We will address these problems and the 
studies of forces and stresses in future publications. 
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VI. APPENDIX 
A. Review of the USPP formulation 



dopotentials can be extremely smooth (i.e. "ultrasoft" ) . 
In order to restore the correctly normalized electron 
charge, the charge density is written as: 



In the USPP approach 34 the constraint of norm- 
conservation on the pseudo wavefunctions is relieved, so 
that the pseudo wavefunctions and the associated pseu- 
I 



KkMI 2 + £ « k |/3 a > Qaa' (r) {fie |Ck) 



(18) 



Here a = [n, I, to, R] and a' = [n', l\ to', R] are sets of orbital quantum numbers and atomic positions R, and (3 a (r) 
are the USPP projector functions. 34 The first term within the square brackets is the ultrasoft charge, and the second 
term the "augmented" charge, that is the portion of the valence charge which is localized within the atomic core radii 
and restores the normalization of the total charge. The atomic charges Q aa > are, by construction, 



Q aa ,(v) = ^(r) - C*(r) ^?», 



lAE 



PS, 



kPSf 



(19) 



where <j)^ E and <p^ s are atomic all-electron and pseudo wavefunctions, respectively. The release of norm-conservation 
leads to the generalized Kohn-Sham equations: 



-V 2 + V LOC (r) + £ \M D «<*> (P*> I + Mr) + V^ c (v) ^ k (r) - < k S< k (r) 



(20) 



r 



where Vloc{ y ) is the local part of the pseudopotential, 
D a aa , is the non-local part and S is the overlap matrix 
which generalizes the orthonormality condition: 



S=l + J2 I/ 3 ") 1- a > (Pa 



(21) 



Here q aa > are the integrals of the augmented charges 
Qaa'(r), and < ip n k\S\ip n 'k' >= ^n,«' h,k'- 

Finally, the non-local pseudopotential projector is 
made up of two contributions: 



u aa' 



= D aa > + 



-J dr (V LOC (r) + V£ xc (r)) Q aQ ,(r). (22) 



The first term on the right side of Eq.22 is the usual 
Kleinman-Bylander projector, and contributes to the 
'bare' pseudopotential (i.e. it is calculated within the 
atomic reference configuration) . The second term is spe- 
cific to the USPP formalism, and represents the action 
which the local and screening potentials exert on the aug- 
mented charges. Since this term depends on V§ xc , it has 
to be updated during the self-consistency cycle. 

In order to ensure better transferability, two atomic 
reference states, corresponding to different energy val- 
ues, are usually included in the projector for each angular 
quantum number. As a consequence, the USPP projector 
contains non-diagonal terms (a, a'), where a = [yi, I, to), 
and a' — (v[ 1 1, to), and v\, v{= 1,2. The atomic reference 
eigenstates do not need to correspond to bound, normal- 
ized solutions of the free atom, but may be unphysical 
eigenstates of the Schr odingcr equation, useful to extend 
the pseudopotential transferability into a larger energy 



range. Thus the atomic pseudo wavefunctions may di- 
verge at large r, but the projector functions /3 Q (r) and 
the matrix D aa > are always short-ranged and well-defined 
by construction. 



B. Pseudo-SIC within USPP 

The USPP implementation of the pscudo-SIC requires 
some generalization of the formalism described in Sec- 
tions III A and IIIB. The charge densities nf of the 
(pseudo) atomic orbitals & are: 

<(r) =v1 (|^(r)| 2 +£(^|/3 Q )Q Qa ,(r)(/3 Q ,|^)) (23) 

\ aa.' / 

and the occupation numbers pi become: 

nk 



(24) 



Furthermore, at variance with the norm-conserving 
case, the non-local part of the USPP depends self- 
consistently on the screening potential itself (see Eq.22). 
As a result it also must be self-interaction corrected. The 
SI part of the non-local USPP is given by: 



1G 



VZs - EE \M JdW^ xc K{v)-l]Q aa ,(v) \ (f3 Q 

1. rvrv 1 \ " / 



(25) 



Thus, the pscudo-SIC KS equations finally are: 



Vloc + Vj 



HXC 



Ei&*>^< (0°'\-(yi IC +vss) 



ICk) = < k 5|Ck), 



(26) 



where Vg IC is given by Eqs.5,12, and 13. 

The total energy is the same as that given in Eq.17, except for the fifth term which now is: 

nk,ir 



(27) 
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